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1.  INTRODUCTION 


The  objective  of  this  work  was  to  develop  adaptive  finite  element  methods  for 
nonlinear  structural  dynamics.  Adaptive  methods  are  particularly  promising  for  nonlinear 
problems  involving  failure,  because  in  failure  and  near-failure  states  of  structures,  three 
phenomena  are  predominant; 

1.  buckling 

2.  shear  banding 

3.  fracture. 

All  of  the  above  phenomena  are  associated  with  localization  of  the  deformation,  by 
which  is  meant  the  development  of  large  strains  in  small  regions  of  the  structure,  which  is 
accompanied  by  large  gradients  in  the  strain.  For  example,  while  strains  are  rather 
distributed  in  elastic  buckling,  once  plasticity  develops  a  large  part  of  the  deformation  of  a 
beam  or  shell  usually  occurs  over  narrow  zone  called  a  hingeline.  Shear  banding  is  a  result 
of  strain  softening  material  behavior  and  is  also  associated  with  very  narrow  x  ands  of 
highly  strained  material.  In  specimens  ranging  in  size  from  0.1 1  to  1.0  meter,  shear  band 
widths  are  of  the  order  of  10  to  100  microns.  In  fracture,  high  strain  gradients  occur  at  the 
crack  tip,  and  in  addition  the  displacement  field  is  discontinuous  behind  the  crack  tip. 

Because  of  the  localization  of  deformation  in  these  problems,  it  is  clear  that  uniform 
meshes  will  be  quite  ineffective.  Although  some  analysts  have  sufficient  intuition  to  guess 
where  the  localization  will  occur,  generally  it  cannot  be  known  a  priori. 

For  these  reasons,  it  seems  that  adaptive  methods  are  essential  in  the  solution  of 
nonlinear  problems.  However,  perusal  of  the  reviews  of  the  field  of  adaptive  finite 
elements  by  Noor  and  Babuska  (Ref.  1)  and  Oden  and  Demkowicz(Ref.  2)  shows  that  the 
majority  of  the  work  has  been  devoted  to  linear  problems. 

In  this  work,  adaptive  methods  are  developed  for  the  nonlinear  dynamics  of  shells 
with  both  geometric  and  material  nonlinearities.  The  localization  phenomenon  which  is  of 
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primary  interest  in  this  class  of  problems  is  hingeline  formation,  but  aspects  of  this  work 
should  be  applicable  to  other  localization  phenomena  in  structural  dynamics.  In  order  to 
simplify  data  management,  an  explicit  time  integration  algorithm  was  chosen.  By  focusing 
the  research  on  explicit  time  integration,  the  research  should  be  relevant  to  a  large  class  of 
widely  used  programs  such  as  DYNA3D. 

There  are  three  types  of  mesh  adaptivity: 

1.  r-adaptive; 

2.  h-adaptive; 

3.  p-adaptive. 

In  the  r-method,  nodes  are  relocated  in  order  to  place  most  of  the  resolution  in  the 
areas  where  they  are  needed.  This  corresponds  in  fact  to  Arbitrary  Lagrangian-Eulerian 
methods;  these  were  studied  in  this  research  and  reported  in  Ref.  3.  However,  r-methods 
were  found  to  be  unsatisfactory  for  treating  localization  phenomena  for  3  reasons: 

1.  in  moving  the  nodes  to  areas  of  high  strain  gradients,  the  elements  become 

severely  distorted,  which  compromises  their  effectiveness; 

2.  it  is  difficult  to  obtain  sufficient  resolution  by  simply  moving  nodes; 

3.  in  shell  problems,  it  is  difficult  to  maintain  fidelity  in  the  modeling  of  the 

shape  of  the  shell  as  the  nodes  move. 

In  p-methods,  resolution  is  increased  where  it  is  needed  by  increasing  the  order  (or 
power)  of  the  interpo’ants.  This  method  was  deemed  inappropriate  for  explicit  dynamics 
methods  for  two  reasons: 

1.  it  is  difficult  to  construct  an  accurate  diagonal  mass  for  higher  order  elements 
and  explicit  methods  rely  for  their  efficiency  on  the  avoidance  of  triangulating 
a  nondiagonal  mass  matrix; 

2.  higher  order  elements  results  in  very  small  critical  time  steps  and  rather  noisy 
solutions  result  when  used  with  explicit  integration. 


In  an  h-method,  elements  are  subdivided  where  greater  resolution  is  needed.  We 
will  call  this  process  fission.  In  addition,  in  dynamic  problems  it  is  possible  to  conserve 
resources  by  fusing  those  elements  which  are  no  longer  actively  deforming.  An  advantage 
of  the  method  is  that  it  is  possible  to  allocate  a  large  number  of  unknowns  in  a  small 
subdomain  without  mesh  distortion  or  introducing  higher  order  elements.  The  h-method 
has  the  disadvantage  that  the  critical  time  step  of  the  mesh  decreases  dramatically  after 
several  levels  of  fission,  but  this  can  be  overcome  by  using  different  time  steps  in  different 
parts  of  mesh. 

Based  on  these  factors,  we  concentrated  on  the  h-method.  This  method  briefly 
summarized  in  Section  2.  Section  3  gives  some  examples  obtained  by  these  methods  and 
Section  4  summarizes  the  research  and  discusses  future  directions. 

2.  SUMMARY  OF  METHODOLOGY 

The  methodology  is  described  in  Refs.  3  and  4.  In  this  Section,  some  of  the  issues 
are  briefly  reviewed.  An  explicit  finite  element  program  for  the  nonlinear  transient  analysis 
of  shells  was  used  as  the  framework  for  the  research.  A  one-point  quadrature  shell  element 
which  Hallquist  in  DYNA3D  manuals  calls  the  Belytschko-Tsay  element(Ref.  5)  was  used. 
We  encountered  some  difficulties  with  this  element  so  they  were  rectified  as  described  in 
Appendix  A.  As  indicated  in  Fig.  1,  the  h-method  consists  of  two  processes: 

1.  fission  of  a  quadrilateral  element  into  4  elements  when  more  refinement  is 
needed; 

2.  fusion  of  four  quadrilaterals  into  a  single  quadrilateral  when  activity  has 
ceased  in  a  subdomain. 

The  fission-fusion  processes  are  triggered  by  error  criteria.  In  this  work  we  have 
studied  two  error  criteria: 

1.  the  error  in  the  flexural  dissipation  rate  due  to  errors  in  the  approximation  of 
the  transverse  displacement; 
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Figure  1, 


Depiction  of  fusion  and  fission  processes  in  h-adaptivity . 
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2.  the  error  in  energy  due  to  one-point  quadrature. 

It  has  also  been  attempted  to  combine  the  two  methods  by  using  the  rates  of  dissipation 
associated  with  the  two  types  of  errors. 

The  studies  have  shown  that  the  first  error  criterion,  error  in  flexural  energy 
dissipation,  as  measured  by  the  energy  dissipated  by  the  discontinuity  in  slope  between  two 
elements,  is  the  more  effective.  Furthermore,  it  has  been  difficult  to  derive  error  criteria 
which  are  effective  for  the  combined  energies. 

Another  aspect  of  adaptive  methods  in  transient  analysis  is  that  error  criteria  should 
only  be  applied  at  selected  intervals  of  the  time  integration  process.  If  the  error  criteria  and 
adaptation  is  attempted  every  time  step,  then  the  fission-fusion  process  becomes 
oscillatory,  i.e.  fission  in  an  element  or  group  of  elements  is  rapidly  followed  by  fusion. 
This  occurs  because  the  fission  process  reduces  the  error  in  a  group  of  elements 
dramatically,  and  if  the  errors  in  these  elements  are  now  compared  to  the  remainder  of 
elements,  the  error  tends  to  fall  near  the  bottom. 

One  technique  which  has  been  successful  in  eliminating  this  churning  is  to  consider 
the  error  criteria  over  groups  of  elements.  In  other  words,  instead  of  using  an  error 
indicator  that  pertains  to  a  single  element,  the  error  is  evaluated  over  a  group  of  elements 
and  then  normalized  by  the  volume  of  elements. 

A  second  technique  we  have  found  very  useful  is  to  limit  the  adaptation  process  to 
selected  time  steps.  For  example,  in  the  problems  that  have  been  studied  here,  limiting  the 
adaptation  process  to  every  50  time  steps  increases  the  efficiency  and  improves  the 
conditioning  of  the  solution.  The  latter  occurs  because  adaptation  tends  to  introduce  noise 
in  the  numerical  solution.  However,  when  the  number  of  time  steps  between  adaptation  is 
moderately  large,  backtracking,  where  the  integration  from  the  previous  adaptation  is 
repeated  with  the  new  mesh,  is  useful. 

One  aspect  of  an  h-adaptive  method  that  requires  considerable  work  is  the  data 
structure.  It  is  necessary  to  keep  track  of  the  elements  from  which  an  element  originated. 
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and  the  siblings(elements  formed  by  the  same  fission).  In  multilevel  adaptivity,  this 
information  must  be  available  for  several  generations.  To  accomplish  this  task,  a  data 
structure  similar  to  that  of  Devloo,  et  al.(Ref.  6)  was  used.  However,  it  was  modified  so 
that  it  is  more  easily  vectorizable. 

It  is  not  clear  whether  h-adaptive  methods  will  be  successful  in  resolving  shear 
bands  in  metallic  structures.  Shear  bands  in  a  specimen  of  size  10" 1  m  are  often  of  the 
order  of  10  microns(10"5  m).  If  it  is  desired  to  capture  any  of  the  structure  of  the  shear 
band,  then  meshes  with  at  least  5  to  10  elements  across  the  band  will  be  needed.  This 
entails  elements  with  h  =  10"6.  Even  if  elements  of  this  size  are  used  only  on  the  band 
itself,  more  than  106  elements  would  be  needed(10  element  across  the  band  by  105 
elements  along  the  band  plus  elements  for  the  rest  of  the  mesh).  Therefore,  some  type  of 
directional  p-refinement  or  methods  such  as  the  spectral  overlay  may  be  more  suitable. 

3.  NUMERICAL  EXAMPLES 

In  this  section,  some  of  the  numerical  results  obtained  by  the  adaptive  methods  will 
be  described.  Since  closed-form  solutions  are  not  available  for  nonlinear  transient 
problems,  two  types  of  comparisons  are  used  for  the  adaptive  solutions: 

1.  numerical  results  obtained  by  finer  meshes; 

2.  experimental  results. 

The  first  example  is  a  cylindrical  panel  which  is  impulsively  loaded  over  the  top  by 
an  explosive  sheet  in  the  area  indicated  in  Fig.  2.  Deformed  mesh  plots  for  the  finer 
adaptive  mesh  are  shown  at  various  times  in  Figs  3  and  4.  Here  the  incremental  transverse 
bending  energy  is  used  for  the  fission-fusion  criterion.  It  can  be  seen  that  after  0.0125 
msec,  the  crown  settler  downward  like  a  plateau  and  the  fissioning  process  migrates 
laterally  towards  the  line  where  the  curvature  is  maximum.  During  this  time,  the  crown 
moves  down  in  a  frozen  plateau-like  state.  After  that,  the  crown  develops  a  convex 
curvature  when  viewed  from  above,  and  the  elements  in  the  crown  are  again  fission-fusion. 
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Figure  3.  Initial  and  deformed  mesh  for  the  cylindrical  panel  with  multi¬ 
level  adaptivity. 
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Figure  4.  Deformed  meshed  for  Che  cylindrical  panel. 
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which  is  a  tendency  that  needs  to  be  fixed:  it  is  probably  due  to  the  fact  that  the  incremental 
work  is  quite  small  in  the  later  stages  because  most  of  the  deformation  has  taken  place,  so 
the  incremental  work  is  quite  uniformly  distributed,  allowing  the  fission-fusion  process  to 
be  triggered  by  small  oscillations  in  the  solution. 

The  time  histories  of  the  displacement  at  two  points  obtained  by  the  h-adaptive 
method  are  compared  to  the  experimental  results  in  Fig.  5  and  6.  Experimental  results  have 
been  obtained  for  this  shell  by  Morino,  Leech,  and  Witmer(Ref.  7),  who  report  a  maximum 
deflection  of  1.27  in.  at  point  A.  The  first  fixed  and  adaptive  meshes  yield  maximum 
deflections  of  1.17  and  1.24  in.,  respectively.  These  results  show  that  even  two-level 
adaptivity  is  quite  successful. 

To  illustrate  the  effectiveness  of  the  adaptive  method  more  clearly.  Table  1  gives  the 
results  for  the  maximum  displacements  at  points  A  and  B  for  various  uniform  meshes  and 
the  adaptive  meshes.  The  finest  meshes  were  run  on  a  Connection  Machine(CM) 
implementation  of  the  program  and  would  not  be  feasible  on  a  CRAY.  It  is  apparent  with 


Mesh  for  half-panel 

Description 

Maximum  Displacement 
(inch)  at 

n  =  -6.2  n  =  -9.42 

6x16 

uniform 

1.001 

0.494 

8x16 

uniform 

1.106 

0.526 

10x20 

uniform 

1.126 

0.527 

12x24 

uniform 

1.146 

0.531 

16x23 

uniform 

1.194 

0.585 

512 

3  level  adaptive 

1.234 

0.606 

64x128 

uniform  CM 

1.264 

0.636 

64x256 

uniform  CM 

1.261 

0.631 

128x256 

uniform  CM 

1.264 

0.636 

Experimental 

B aimer  and  Winner  (1964) 

1.28 

0.68 

Table  1.  Cylindrical  Panel  Results 
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these  extremely  fine  meshes,  a  converged  solution  has  been  obtained.  The  converge 
solution  does  not  agree  exactly  with  the  experimental  solution,  and  while  the  ultimate  test  of 
a  computation  is  agreement  with  experimental  results,  in  nonlinear  problems  it  is  difficult  to 
distinguish  errors  arising  from  deficiencies  in  the  constitutive  models,  boundary  modeling 
and  loading  from  errors  due  to  insufficient  mesh  resolution.  In  the  CM  runs,  resolution  is 
apparently  no  longer  an  issue,  so  these  results  can  be  considered  a  benchmark  for  this 
problem.  It  is  interesting  to  observe  how  slowly  uniform  meshes  converge  to  this 
benchmark.  Even  a  12  x  24  mesh(for  half  of  the  panel)  is  more  than  10%  in  error;  many 
results  reported  in  the  structure  involve  even  coarser  meshes.  On  the  other  hand  an 
adaptive  mesh  with  the  same  member  of  elements  is  only  4%  in  error. 

The  second  example  is  a  hollow,  cylindrical  column  which  is  subjected  to  a 
compressive  axial  load.  This  problem  is  of  interest  because  it  exhibits  both  global  and  local 
buckling,  the  latter  resulting  in  buckling  of  the  cross  section.  Numerical  results  and 
experimental  results  are  reported  for  this  problem  by  Kennedy,  Belytschko,  and  Liu(Ref. 
8). 

The  cylinder  is  loaded  by  prescribing  an  upward  velocity  of  500  in/sec  to  the 
bottom  nodes  of  the  model,  with  the  top  fixed.  To  trigger  the  lateral  buckling  mode,  an 

imperfection  given  by 

-  .  271Z 

Ax  =  0.01  sin  — j~ 

where  fis  the  length  of  the  column  and  z  is  the  coordinate  along  the  axis  of  the  column,  is 
added  to  the  x-coordinate  of  all  nodes. 

The  problem  is  solved  with  multi-level  adaptivity  and  the  transverse  energy 
criterion.  The  pattern  of  adaptivity  is  shown  in  Fig.  7.  Initially,  the  fission  process  is  only 
one  level  and  occurs  at  the  eventual  extreme  of  the  buckling  wave.  The  important  point  to 
note  is  that  the  fission  process  then  focuses  at  the  points  of  local  buckling. 
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Figure  7. 


Deformed  adaptive  meshes  for  the  cylindrical  column  with  multi¬ 
level  adaptivity. 
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Another  class  of  example  we  have  studied  is  shear  banding.  H-adaptivity  was  not 
successful  here  so  we  developed  a  new  method  which  is  in  a  sense  p-adaptive.  In  this 
method,  the  resolution  is  increased  by  overlaying  the  finite  element  mesh  with  a  spectral 
interpolant  in  a  zone  that  coincides  with  the  shear  band.  The  results  for  a  tensile  specimen, 
shown  in  Fig.  8,  are  given  in  Figs  9  to  11.  The  noteworthy  results  are  the  complex 
structure  of  the  shear  band.  The  sequential  times  A,  B  and  C  are  indicated  in  these  figures. 
Most  of  the  evolution  of  the  shear  band  takes  place  between  B  and  C.  Note  that  the  shear 
band  has  considerable  structure,  and  as  can  be  seen  from  Fig.  1 1,  its  width  is  of  the  order 
of  several  hundred  microns.  This  is  larger  than  observed  physically  but  much  smaller  than 
the  element  size.  We  have  since  learned  that  the  width  is  determined  by  the  scale  of 
imperfections. 

4.  SUMMARY  AND  DISCUSSION 

The  results  have  shown  that  an  adaptive  mesh  can  dramatically  improve  accuracy 
for  a  transient,  nonlinear  finite  element  solution.  Comparison  of  results  for  the  cylindrical 
panel  with  a  converged  solution  obtained  on  a  massively  parallel  computer  on  a  very  large 
mesh  of  16K  elements  show  that  the  h-adaptive  method  can  obtain  similar  accuracy  with 
528  elements,  whereas  a  uniform  mesh  of  this  size  is  in  error  in  the  maximum  displacement 
by  almost  10%. 

A  key  finding  of  this  work  is  that  the  error  criteria  generally  used  in  linear  analysis 
are  not  well  suited  for  problem  of  localization  in  nonlinear  mechanics.  The  difficulty  arises 
because  linear  error  criteria,  such  as  that  of  Zienkiewicz  and  Zhu  consider  the  difference 
between  a  C°  least-square  fit  to  the  stresses  and  the  C’1  stress  field  obtained  by  the  finite 
element  scheme;  the  error  criteria  of  Babuska  relate  the  error  to  the  jumps  in  the  stress  fields 
at  element  interfaces,  but  give  similar  results.  In  localization,  such  as  buckling  and  shear 
banding,  the  stress  field  is  very  homogeneous,  so  the  difference  between  the  C°  and  C’1 
fields(and  the  jumps  in  the  stresses  at  element  interfaces)  are  small.  Thus  these  error 
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Figure  10.  The  maximum  effective  plastic  strain  in  the  shear  band  as  a 
function  of  the  displacement  of  the  ends  of  the  specimen. 


Figure  11. 


Effective  plastic  strain  in  a  cross-section  across  the  shear 


criteria  will  not  be  successful  in  determining  where  to  refine  the  mesh.  It  appears  the 
development  of  refinement  criteria  for  such  problems  may  be  accomplished  most  effectively 
if  the  physics  of  the  phenomena  which  are  expected  are  used  as  guideline.  Thus  in  this 
study,  where  hingeline  formation  was  of  interest,  we  used  the  flexural  energy  as  an  error 
indicator.  Although  such  approaches  would  not  be  suitable  for  completely  automatic 
adaptive  computation  of  nonlinear  problems,  it  has  appeal  for  skilled  analysts. 

In  this  study,  the  localization  phenomena  considered  was  the  formation  of  elastic- 
plastic  hingelines  associated  with  buckling.  The  application  of  these  techniques  to  other 
localization  phenomena  requires  further  work.  In  particular,  error  criteria  appropriate  for 
other  types  of  localization  need  to  be  devised.  The  error  criteria  developed  here  are  not 
universally  applicable  to  nonlinear  problems.  Considerable  additional  work  is  needed  on 
error  criteria  if  adaptive  methods  are  to  be  usefully  applied  to  nonlinear  structural  problems 
with  localization.  Since  accurate  resolution  of  localization  is  the  major  reason  for  applying 
adaptive  methods  to  nonlinear  structural  problems,  unless  error  criteria  are  developed, 
adaptivity  will  not  be  successful. 

To  avoid  excessive  churning  of  the  fission-fusion  process,  time  delays  had  to  be 
incorporated  in  the  judgment  process.  Thus,  fission  and  fusion  are  executed  only  when 
indicated  by  two  or  more  consecutive  judgment.  Nevertheless,  churning  becomes  a 
problem  with  the  incremental  dissipation  criterion  in  the  later  stages  of  impulsively  loaded 
problems  when  the  work  on  the  system  decreases;  rational  methods  of  controlling  it  need  to 
be  developed. 

The  h-adaptive  procedure  is  limited  in  its  ability  to  resolve  the  subdomains  of 
maximum  deformation  by  the  fact  that  the  parent  element  configuration  is  fixed.  Therefore, 
hinge  lines  which  occur  at  angles  relative  procedure  to  the  mesh  lines  may  not  be  captured 
effectively.  However  the  h-method  appears  to  be  the  best  compromise  between  simplicity 
and  effectiveness  in  the  solution  of  nonlinear  structures  by  explicit  methods.  The  results 
we  have  obtained  show  that  these  adaptive  schemes  are  capable  of  achieving  substantial 
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improvements  in  accuracy  for  a  given  computational  effort.  Generally,  an  adaptive  mesh  is 
capable  of  achieving  the  same  accuracy  with  less  than  half  of  the  computational  resources. 
The  fission  process  tends  to  take  place  in  the  subdomains  where  the  maximum  deformation 
occurs. 

In  this  study,  an  h-method  was  selected  because  of  its  advantages  in  an  explicit 
integration  algorithm.  However,  further  study  on  quadratic  element  technology  has 
indicated  that  effective  9  node  elements  can  be  developed  with  4  quadrature  points.  It  has 
also  been  found  that  by  iterating  on  the  off-diagonal  terms  of  the  mass  matrix  only  once  or 
twice,  it  is  possible  to  obtain  good  solutions  to  transient  problems  with  these  higher  order 
elements  without  triangulating  the  complete  consistent  mass.  Therefore  p-method  deserves 
further  consideration  in  these  problems. 
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APPENDIX  A 


I  INTRODUCTION 

Four-node  quadrilateral  shell  elements  with  one  quadrature  point  in  the  mid-surface  have 
become  widely  used  in  programs  with  explicit  dme  integration.  The  first  of  these  elements  was 
described  by  Belytschko  and  Tsay  (Ref.  3);  also  see  Belytschko,  Lin  and  Tsay  (Ref.  4).  This 
element  is  used  in  DYNA3D,  PAMCRASH  and  other  programs  developed  for  crashworthiness 
studies.  Hallquist  (Ref.  7)  also  adapted  the  Hughes  and  Liu  (Refs.  8  and  9)  shell  element  to 
these  programs  by  adding  an  hourglass  control  similar  to  that  in  Ref.  4. 

The  major  objective  in  the  development  of  the  Belytschko-Tsay  element  was  to  attain  a 
convergent,  stable  element  with  the  minimum  number  of  computations.  For  this  reason,  the 
element  uses  bilinear  isoparametrics  with  one  quadrature  point  in  the  mid-plane  when  the 
material  is  elastic.  If  the  material  is  nonlinear,  several  quadrature  points  are  used  through  the 
thickness  at  a  single  mid-plane  point.  Because  this  element  with  one-point  quadrature  would  be 
rank  deficient,  i.e.,  it  would  possess  so-called  "hourglass  modes"  or  "spurious  singular  modes", 
an  hourglass  control  is  added.  This  hourglass  control  is  orthogonal  to  all  linear  fields  (see 
Belytschko  and  Tsay  (Ref.  3)  and  Belytschko,  Lin  and  Tsay  (Ref.  4)),  so  the  consistency  for 
linear  fields  is  not  impaired. 

Because  of  the  emphasis  on  speed,  several  shortcuts  were  made  in  formulating  the  element 
equations.  On  the  whole,  the  element  has  performed  quite  well,  but  it  has  two  shortcomings: 

1.  It  performs  poorly  when  warped,  and  in  particular,  it  does  not  correctly  solve  the 
twisted  beam  problem. 

2.  It  does  not  pass  the  quadratic  Kirchhoff-type  patch  test  in  the  thin  plate  limit. 

The  latter  shortcoming  is  shared  by  the  Hughes-Liu  element,  and  its  importance  was  not 
realized  until  recently. 

In  this  paper,  modifications  to  the  Belytschko-Tsay  element  which  overcome  these 
drawbacks  are  described.  The  first  shortcoming  is  eliminated  by  adding  terms  into  the  strain- 
displacement  equations  which  couple  the  curvatures  to  the  translations:  these  terms  are  shown 
to  be  essential  to  obtaining  the  correct  response  when  the  shell  is  twisted  and  the  mid-plane 
shape  depends  strongly  on  the  bilinear  term. 

The  second  shortcoming  is  eliminated  by  adding  a  nodal  projection  to  the  shear  calculation. 
The  concept  of  a  projection  operator  to  improve  element  performance  was  originated  by 
Belytschko,  Stolarski  and  Carpenter  (Ref.  5)  and  studied  in  a  quadrilateral  by  Stolarski, 
Carpenter  and  Belytschko  (Ref.  15).  The  particular  projection  used  here  is  similar  to  that 
developed  independently  by  Hughes  and  Tezduyar  (Ref.  10)  and  MacNeal  (Ref.  12). 

In  addition,  we  describe  an  implementation  of  the  patch  test  for  explicit  dynamic  programs. 
The  patch  test  is  usually  defined  in  terms  of  a  static  analysis,  which  is  not  possible  with  a 
program  that  includes  only  explicit  time  integration.  The  patch  test  described  here  can  be  used 
directly  in  explicit  programs  without  any  modification,  so  it  should  prove  useful. 

This  element  formulation  is  based  on  the  "resultant  stress  theory"  (Liu,  et  al.  (Ref.  11)  and 
Stanley  (Ref.  14)).  The  starting  point  is  the  degenerated  continuum  (DC)  approach  to  shells 
which  describes  the  shape  and  kinematics  (Ref.  8).  The  strains  are  then  expressed  in  terms  of 
membrane  strains  and  curvatures.  The  advantage  of  resultant  stress  theories  over  DC  shell 
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elements  is  that  the  number  of  computations  is  substantially  smaller,  which  provides  significant 
speed  advantages  in  explicit  computer  programs. 

The  paper  is  organized  as  follows:  In  Section  2,  the  starting  point,  the  DC  shell  geometry 
and  interpolation  are  reviewed  along  with  the  corotational  coordinate  system  used  in  this 
element.  Section  3  describes  the  new  kinematic  relations;  two  methods  are  presented,  one 
requires  a  knowledge  of  the  pseudonormal  vectors,  the  second  does  not.  Section  4  describes  the 
shear  projection.  Section  5  describes  the  implement?  .ion  of  the  element,  giving  the  rate-of- 
deformation-velocity  equation.  Section  6  gives  some  numerical  results  which  contrast  the 
performance  of  this  element  with  the  earlier  element;  a  recipe  for  performing  the  patch  test  in 
explicit  programs  is  also  given. 


2.  GEOMETRY  AND  INTERPOLATION 


In  DC  shell  theories,  the  coordinates  of  the  shell  are  given  by: 

me 

x=4X  kop(1+0  +  x^a-oH^r!) 
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Thus,  if  we  define  the  pseu  ionormal  at  each  node  I  by: 


(la) 

(lb) 

(lc) 

(2) 


where  h{  is  the  distance  between  the  two  nodes  (or  a  pseudo-thickness  at  node  I),  and  define  the 
nodes  on  the  mid-surface  by: 


then: 


where: 


and 


=  i-(x|°P+x!>°1) 

(3) 

me 

=  Xx-  +  Cp 

(4a) 

1=1 

me 

=  X  (xrKPi)Ni(£,Ti) 

(4b) 

i=i 

«rt| 

II 

*\9 

(5) 

-Z?.- 


(6) 


X">=]T  X^.T]) 


1=1 


In  the  above,  me  is  the  number  of  element  nodes,  Ni  are  the  shape  functions,  and  h  is  the 
pseudo-thickness.  Uppercase  indices  refer  to  nodal  numbers.  The  superscripts  "top"  and  "bot" 
refer  to  the  top  and  bottom  nodes  which  describe  the  actual  continuum,  whereas  the  superscript 
"m"  refers  to  the  mid-surface;  nodal  coordinates  without  superscripts  pertain  to  the  midsurface. 
The  geometry  is  shown  in  Fig.  1. 

The  element  corotational  system  (x,y,z)  is  constructed  as  shown  in  Fig.  2.  The  mid-points 
of  the  sides  are  connected  by  lines,  rac  and  rbd>  as  shown.  The  direction  of  the  z  coordinate, 
which  corresponds  to  the  unit  vector  e3,  is  then  obtained  by: 


e,  _  racxrbd 
1 1  racxrbdl  I 


(7) 


The  precise  orientation  of  the  two  other  unit  vectors  is  usually  not  important;  we  choose: 

ei  =  rac/||rac||  (8) 

e2  =  e3xei  (9) 


The  choice  of  the  unit  normal  e3  is  also  not  critical,  but  this  particular  choice  has  an 
interesting  consequence.  Because  the  two  lines  rac  and  ije  in  the  surface  of  the 
isoparametric  element,  the  vector  e3  is  exactly  normal  to  the  mid-surface  at  the  origin  of  the 

reference  plane.  By  contrast,  when  e3  is  the  cross-product  of  the  diagonals  of  the  element  as  in 

Belytschko,  Lin  and  Tsay  (Ref.  4),  then  e3  is  not  exactly  normal  to  the  element  surface,  which 
complicates  some  of  the  subsequent  developments.  The  practical  effects  are  probably  not 
important. 

The  mid-surface  of  the  shell  is  given  by: 

4 

(io) 
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where  Ni(x,h)  are  the  usual  bilinear  isoparametric  shape  functions,  which  are  given  by: 


Ni  =  7  (!+££)(  1+T)ir|) 
4 

The  velocity  of  the  shell  is  obtained  from  Eq.  4: 


(ID 
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Superposed  dots  throughout  this  paper  indicate  material  time  derivatives. 
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3.  KINEMATICS 


3.1  General  Description 

The  rate-of-deformation  (or  stretching)  tensor  in  the  corotational  system  is  given  by: 
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(3vi 
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To  evaluate  this  tensor,  we  need  to  obtain  the  derivatives  of  the  velocity  field.  From  Eqs.  4  and 
12  it  follows  that  the  velocity  field  is  given  by: 
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Using  implicit  differentiation,  the  derivatives  of  the  shape  functions  are  given  by: 
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where  J  is  the  Jacobian.  Now  using  Eq.  4  to  express  the  derivatives  of  x,y  with  respect  to  (x,h), 
and  noting  that  Eq.  4  holds  in  any  coordinate  system,  we  obtain: 
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As  described  in  Belytschko,  Wong  and  Stolarski  (Ref.  6),  the  terms  in  J  which  are  linear  or 
higher  order  in  C,  in  Eq.  17c  have  little  effect  on  element  performance,  so  the  only  term  involving  £ 
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which  will  be  retained  is  the  second  term  in  Eq.  16.  Thus,  using  Eqs.  13,  14  and  16  gives  the 
following  stretching-velocity  relations: 
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At  the  quadrature  point,  x  =  h  =  0,  the  matrix  bn  is  that  given  in  Belytschko,  Lin  and  Tsay 
(Ref.  4): 
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Two  methods  have  been  developed  to  evaluate  bu,  the  terms  which  couple  curvatures  to 
translations.  The  two  methods  are  distinguished  as  follows: 

1.  Method  p  involves  the  direct  evaluation  of  the  second  term  in  Eq.  16;  it  requires  the 
p  vectors  to  be  available  at  all  nodes,  which  is  not  the  case  for  the  Belytschko-Tsay 
element. 

2.  Method  z  involves  an  estimate  of  these  terms  based  on  the  assumption  of  normality 
of  p  to  the  surface  of  the  shell;  the  p  vectors  need  not  be  available. 

3.2  Method  p 

In  Method  p,  the  terms  such  as  Px.n-are  evaluated  directly.  The  computation  is  quite  simple; 
from  Eq.  4: 


4  4 

Px.ti  =  X  PSjNj.n  =  \  X  (2°) 

j=i  4  J=i 

where  the  last  step  follows  from  differentiation  of  Eq.  11.  Hence,  using  a  similar  formula  for  Px£> 
it  follows  that: 
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|bxl|  1  ^  /  ^iPyKHK^HlPyK^sK  | 

Wi j  16J  1^1  \  -SlTCldlK+TllP^K  / 


Py2“Py4  Py3~Pyt  Py4"Py2  Py 1 ~Py 3 
PJcPfc  Pxl'K3  Pfc'Pfc  K3-Kl 


1=1  1=2  1=3  1=4 


(21a) 

(21b) 


3.3  Method  7. 

In  Method  z,  the  normal  to  the  mid-surface  is  first  determined.  The  derivatives  of  this  vector 
then  give  b^.  To  obtain  the  normal,  we  start  with  an  expression  for  the  mid-surface  based  on  the 
bilinear,  isoparametric  field  given  by  Belytschko  and  Bachrach  (Ref.  2): 


4 

z™  =  X  (si+xbxi+ybyi+^yOzi 
1=1 


si  =[1,1,  1,1] 


X  hJ*J 

J=l 


by] 


hl  =  [+l,-l,+l,-l] 


(22) 

(23a) 

(23b) 

(23c) 

(23d) 


where  bii  are  defined  in  Eq.  19. 

The  normal  to  the  surface,  p,  is  obtained  by  taking  the  gradient  to  the  surface  described  by 
Eq.  22,  which  gives: 


=  --Ly 


zi 


1=1 


I  bxi+(£n)  $Yi  \ 

j  byi+(^T|),9Yi  j 


(24a) 


where: 


p**(l+z|+z$)1/2  (24b) 

At  the  origin  (^'H).x  -  (^'H).y  =  0,  because: 

(^ri),x  =  ^rj.x  +  ^  =  0  (25) 
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As  described  earlier,  p  by  construction  is  normal  to  the  x-y  plane  at  the  origin,  so  from  Eqs.  24a 
and  25,  it  follows  that: 

4  4 

X  bxizi  *  X  hyizi  =  0  (26) 

1=1  i=i 

Therefore  p*  =  1  at  the  origin  of  the  reference  plane,  i.e.  at  the  quadrature  point. 

Taking  the  derivatives  of  pj  and  pp  with  respect  to  x  and  h  (and  neglecting  the  terms  related 

to  p  ^  and  p,^,  which  can  be  shown  to  be  small)  gives: 


re.$  =  =  -vu 

Px,tj  =  "z,xn  =  -Zy4.x 

Py,£  =  =  -ZyT1.9 

Py,r|  =  ‘z.yn  =  ‘^.y 


(27a) 

(27b) 

(27c) 

(27d) 


where: 


V-''  /V 

Zy  =  yTZ  =  2,  YIZI 


(27e) 


Since: 


where: 


$.*  iy  j  =  lf  ^  'Xll 

Tl.x  ”H  ,y  J  -y^  X£  J  -y^  x^ 


4yi]  =  T)~y  =  2*  'Hiyr  etc- 


it  follows  from  Eqs.  16,  27  and  28  that: 


kl|  "  16J2  (  l^.y.  j 

2Zyf  X13  X42  x31  x24 


7  /-S  /N  ^ 

a  L  yi3  y42  y3i  y24 


(30a) 


(30b) 


Thus,  this  matrix  involves  the  same  terms  as  the  b  matrix  given  in  Eq.  19. 
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Remark.  Method  z  couples  curvatures  to  translations  only  for  warped  elements,  i.e.  when  the 
nodes  are  not  coplanar  so  that  zg  *  0.  Method  p,  on  the  other  hand,  also  introduces  a  coupling 
for  surfaces  of  single  curvature.  In  the  latter  case,  however,  the  coupling  appears  to  be 
insignificant. 

4.  SHEAR  PROJECTION 

The  shear  strains  are  calculated  by  a  nodal  projection  based  on: 

e1n=i-(eii+eU+-L-(v2j-vZI)  (3i) 

2  Lu 

where  the  superscript  I  refers  to  side  I  and  the  subscript  n  refers  to  a  component  normal  to  side  I; 
see  Fig.  3. 

The  transverse  shears  then  are  given  by: 

Yxz  =  £  Ni(S  ,  ri)e?,  (32a) 

i=i 

Yyz  *  -  X  -  *1  )0a  (32b) 

1=1 

The  transverse  shears  do  not  depend  on  w  because  after  the  projection,  w  is  considered  to  have 

vanished  (see  Refs.  5  and  15).  The  terms  9a  are  obtained  from  0„  by  the  standard 
transformation: 


®xi  “(eA-ejjtS  +(e£ej)6n 

(33a) 

0yi=(e^e?)e1n+(enK-e?jenK 

(33b) 

where  e-  and  en  are  unit  vectors  defined  in  Fig.  3. 

Evaluating  the  resulting  forms  for  the  transverse  shear  at  the  quadrature  point,  x  =  h  =  0, 
gives: 


I  3-1 
13.1 


b'j 

b;. 


Vzl 

XV 

9xi 

Oyl 


2(xji-xik) 

2(yji-yiK) 


(xjiyji-xixyiK)  -  (xjixji-xkxik) 

(yjiyji-yncyiic)  *  (xjiyji-xiKyiK) . 


—  xv  FT  _  xv  _  Jl 

xji  =  xn/L  yn  =  yn/L 


(34) 


(35) 

(36a) 
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(36b) 


Ln 


*2 

= xft+yft 


5.  IMPLEMENTATION 

The  implementation  of  this  element  closely  parallels  that  described  in  Ref.  4.  One-point 
quadrature  is  used  in  conjunction  with  hourglass  control. 

The  velocity  strains  for  Method  z  in  the  mid-surface  are  given  by  combining  Eqs.  18,  19  and 
30,  which  gives: 


ST  =  (y24Vx13+yi3Vx42)/(2A) 

(37a) 

d^  =  (X42Vyl3+Xi3Vy24)/(2A) 

(37b) 

2dxy  =  (X42Vx13+Xl3Vx24+y24Vyl3+y3lVy24)/(2A) 

(37c) 

The  curvatures  are  given  by: 

/N 

Kx  =  (y249yl3+yi30y42)/(2A) 

(37d) 

+  2  2^Xi3Vxi3+X42Vx24)/A2 

Ky  =  -(X420xl3+Xi39x24)/(2A) 

(37e) 

+  2  z7(y13Vy13+y42Vy24)/A2 

✓S  A  /\  />» 

2lCxy  =  (X420yi3+Xi30y24-y240xl3-y310x24)/(2A) 

(370 

+  2  Zy(xi3Vyl3+X42Vy24+yi3vxl3+y42vx24 

The  total  velocity  strains  are  then  computed  by: 

dx  =  dx  + 

(38) 

dy  =  dy  +  ZXy 

(39) 

'T  _  2Tn  ~ 
dxy  —  Uxy  +  Z*xy 

(40) 

The  hourglass  strain  rates  are  computed  as  in  Ref.  4;  some  modifications  are  needed  to 
exactly  satisfy  the  patch  test.  The  transverse  shear  velocity  strains  are  computed  as  described 

in  the  previous  section.  The  stresses  sij  and  the  hourglass  stresses  Qi*,  Q^,  Q?,  Qf  and  Qf 
are  then  computed  by  the  constitutive  equation.  The  nodal  force  expressions  then  emanate  from 
the  transpose  of  the  kinematic  relations. 
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If  the  corotational  coordinate  system  xi,X2  is  updated  according  to  the  spin  as  described  in 
Ref.  4,  the  rate  of  the  stress  corresponds  to  the  Green-Naghdi  rate.  The  formulation  thus 
requires  a  constitutive  law  which  relates  the  Green-Naghdi  rate  to  the  corotational  stretching 
tensor  (Eq.  13).  Under  these  conditions,  the  formulation  is  valid  for  large  membrane  strains. 

6.  NUMERICAL  STUDIES 

6.1  The  Patch  Test  for  Explicit  Programs 

The  patch  test  is  of  great  value  in  verifying  the  theoretical  validity  of  elements  and  their 
implementation.  Unfortunately,  it  has  not  been  used  with  explicit  programs  because  the  patch 
test  usually  is  stated  as  a  linear,  static  problem.  This  type  of  problem  is  not  easily  treated  by 
nonlinear,  explicit  transient  programs.  This  section  describes  a  procedure  for  implementing  the 
patch  test  in  an  explicit  program  which  requires  no  modifications  of  most  programs. 

Before  describing  this  procedure,  we  will  describe  the  standard  static  patch  test,  which  will 
clarify  our  implementation  of  the  dynamic  patch  test.  In  the  static  patch  test,  an  irregular  mesh 
such  as  shown  in  Fig.  4  is  considered.  At  the  boundary  nodes,  a  displacement  field  which  is 
consistent  with  a  state  of  constant  strain  is  prescribed.  The  displacements  of  the  interior  nodes 
should  then  be  consistent  with  the  linear  displacement  field  associated  with  this  constant  strain. 

In  the  plane  patch  test,  the  displacements  around  the  periphery  of  the  mesh  are  prescribed 
by: 


ux  =  ctix  +  ot2y  +  0C3  (41a) 

uy  =  0C4X  +  ct5y  +  ct6  (41  b) 

where  btj  are  arbitrary  parameters  set  by  the  user.  Satisfaction  of  the  patch  test  then  requires 
that  in  the  static  solution  for  this  mesh,  the  displacements  of  the  interior  nodes  match  exactly 
those  given  by  Eq.  41. 

In  explicit  programs,  the  implementation  of  a  patch  test  is  hampered  by  the  inability  to  obtain 
an  exact  static  solution.  Moreover,  most  explicit  programs  use  rates-of-deformation  as  a 
measure  of  deformation,  so  prescribing  initial  displacements  will  not  work.  However,  these 
difficulties  can  be  circumvented  as  follows: 

1.  Prescribe  initial  velocities  corresponding  to  Eq.  41  at  all  nodes  of  the  mesh,  i.e. 

ux  =  ai  +  ®2X  +  a3y  (42) 

iiy  =  CC4  +  CX5X  +  ct6y 

2.  Integrate  one  time  step  with  zero  body  forces. 

3.  The  accelerations  should  vanish  at  all  interior  nodes. 

Remark.  The  constants  ai  should  be  small  enough  so  that  any  geometric  non  linearities  brought 
into  play  during  the  time  step  are  insignificant. 

The  concept  underlying  this  test  is  that  if  Eq.  42  is  used  as  an  initial  condition  on  the 
velocities,  then  after  the  first  time  step,  the  mesh  will  displace  into  the  configuration  of  Eq.  41 
with  a;  =  GtjAt.  The  elements  should  then  all  have  the  same  constant  state  of  strain  and  stress, 
so  the  accelerations  at  the  interior  nodes  should  vanish.  At  the  boundary  nodes,  the 
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accelerations  will  not  vanish  because  an  external  force  is  needed  to  generate  this  state  of 
constant  stress. 

For  plate  bending,  the  explicit  patch  test  consists  of  prescribing  iiz,  9X  and  0y  as  follows: 


uz  =  cti  +  ot2X  +  ot3y  +  (X4X2  +  ot5xy  +  a^y2 

(43a) 

9x  =  —  =  oc3  +  oc5x  +  2a6y 

9y 

(43b) 

0y  =  =  a.2  +  ct5y  +  20C4X 

dx 

(43c) 

The  system  is  integrated  one  time  step  with  zero  external  loads.  The  accelerations  of  interior 
nodes  should  then  be  almost  exacly  zero,  with  any  deviations  ascribable  to  geometric 
nonlinearities.  If  the  program  is  a  linear  program,  they  should  be  exactly  zero. 

6.2  Twisted  Beam 

The  twisted  beam  problem  is  described  in  Fig.  5.  In  this  problem,  we  used  5  degrees  of 
freedom  per  node;  6  degrees  of  freedom  still  causes  errors.  The  time  history  of  the  displacement 
at  the  tip  for  the  in-plane  load  is  given  in  Fig.  6  and  is  compared  to  the  Belytschko-Tsay  element. 
As  can  be  seen,  the  latter  diverges  immediately;  without  the  additional  terms  described  here,  the 
beam  has  almost  no  stiffness.  The  present  element  performs  very  well.  Results  have  been 
compared  to  those  obtained  by  the  Hughes-Liu  element;  they  differ  by  less  than  1  percent. 

6.3  Hemispherical  Shell 

The  mesh  for  the  hemispherical  shell  is  shown  in  Fig.  7.  This  mesh  differs  from  the  mesh  in 
the  MacNeal-Harder  test  set  (Ref.  13)  in  that  there  is  no  hole  in  the  top.  Therefore,  while  in  the 
MacNeal-Harder  problem  the  nodes  of  each  element  are  coplanar,  in  this  mesh  they  are  not.  This 
has  a  significant  effect  on  the  performance  of  the  element. 

As  can  be  seen  from  Fig.  8,  the  displacement  time  history  of  the  Belytschko-Tsay  element  is 
quite  erratic  in  this  dynamic  problem.  The  improved  element,  however,  behaves  perfectly. 

6.4  Cylindrical  Panel  Problem 

The  last  test  problem  considered  here  is  an  impulsively  loaded  cylindrical  panel.  The 
problem  description  and  the  mesh  are  given  in  Fig.  9.  Five  integration  points  were  used  through 
the  thickness.  The  number  of  quadrature  points  has  a  sizable  effect  on  the  results;  marked 
stiffening  is  observed  when  increasing  the  number  of  quadrature  points  from  3  to  5. 

A  time  history  of  the  midpoint  deflection  is  shown  in  Fig.  10,  where  it  is  compared  with  the 
experimental  results  of  Balmer  and  Witmer  (Ref.  1)  and  with  results  for  the  Belytschko-Tsay 
element.  The  results  of  the  two  elements  differ  little  because  most  of  the  elements  deform  into 
configurations  in  which  the  nodes  are  coplanar.  It  is  only  when  a  significant  fraction  of  the 
elements  is  warped  that  these  modifications  are  important. 

7.  CONCLUSIONS 
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Additional  terms  have  been  added  to  the  strain-displacement  equations  of  the  Belytschko- 
Tsay  element,  and  the  transverse  shears  have  been  modified  by  a  projection.  These 
modifications  involve  hardly  any  additional  computation  time;  at  mosi,  a  10  percent  increase  in 
the  cost  of  the  element  has  been  observed. 

These  terms  give  dramatic  improvements  in  the  performance  of  the  element  in  the  twisted 
beam  problem  and  for  a  particular  meshing  of  the  hemispherical  shell  problem.  In  other  problems, 
the  differences  are  insignificant.  The  terms  are  important  only  when  the  nodes  of  an  element  are 
not  coplanar. 

An  implementation  of  the  patch  test  for  explicit  programs  has  also  been  described.  This 
element  passes  the  patch  test,  whereas  the  Belytschko-Tsay  and  Hughes-Liu  elements  do  not 
pass  the  patch  test.  Failure  to  pass  the  patch  test  has  significant  ramifications  on  the 
convergence  of  elements  and  their  performance  in  distorted  meshes. 
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O  MIDSIDE  NODES 


Fig  2  Orientation  of  Corotational  Coordinates 
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Fig.  4  Mesh  for  Patch  Test. 
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Fig.  5  Local  Coordinates  for  Rotation  Projection 


Young's  modulus  E  =  29x10 
Poission's  ratio  v  =0.22 


and  Problem  Description  for  Twisted  Beam  Problem 


Fig  7  Displacement  at  Tip  Normalized  to  Exact  Static  Solution 
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radius  =10.0 
thickness  =  0.04 
E  =  6.825x10  7 
v=  0.3 
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Fig.  8  Hemispherical  Shell  Mesh  and  Problem  Description. 


Fig.  9  Displacement  Under  Load  Normalized  to  Exact  Static  Solution. 


thickness  =  0.25 
E  =  4.32x10 
v=  0.33 

weight  =  90.0  per  unit  area 


Fig.  10  Impulsively  Loaded  Cylindrical  Panel 


deflection  (  in.) 


time  (  msec.) 

Fig.  1 1  Vertical  Displacement  of  Cylindrical  Panel  9.42  in.  From  Near  Support 
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